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Abstract. The non-ergodic behavior of the deterministic Fixed Energy Sandpile (DFES), with Bak-Tang- 
Wiesenfeld (BTW) rule, is explained by the complete characterization of a class of dynamical invariants (or 
toppling invariants). The link between such constants of motion and the discrete Laplacians properties on 
graphs is algebraically and numerically clarified. In particular, it is possible to build up an explicit algorithm 
determining the complete set of independent toppling invariants. The partition of the configuration space 
into dynamically invariant sets, and the further refinement of such a partition into basins of attraction for 
orbits, are also studied. The total number of invariant sets equals the graphs complexity. In the case of two 
dimensional lattices, it is possible to estimate a very regular exponential growth of this number vs. the size. 
Looking at other features, the toppling invariants exhibit a highly irregular behavior. The usual constraint 
on the energy positiveness introduces a transition in the frozen phase. In correspondence to this transition, 
a dynamical crossover related to the halting times is observed. The analysis of the configuration space 
shows that the DFES has a different structure with respect to dissipative BTW and stochastic sandpiles 
models, supporting the conjecture that it lies in a distinct class of universality. 
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1 Introduction 

Sandpile models have been introduced in statistical me- 
chanics as prototypes for the Self-Organized Criticality 
(SOC), a concept that has been widely used in order to 
explain the appearance of power law correlations in non- 
equilibrium steady states of self-organizing systems with 
many degrees of freedom |1I2| . As conceived by Bak, Tang, 
and Wiesenfeld (BTW) in their seminal work [3], dissipa- 
tion and external input of grains should be necessary con- 
ditions for the existence of a self-organized critical state, 
j On the other hand, the common numerical technique used 
to study the critical behavior of these models is the finite- 
size scaling, that requires the analysis of larger and larger 
systems [3]. Increasing the system's size implies that, at 
the criticality, larger and larger avalanches are produced, 
i.e. the system's self-sustained activity persists in time. 
Since new grains are added only when the dynamics even- 
tually stops (at least in the original BTW model and 
dissipation is localized at the boundaries, it has been con- 
jectured that the correct infinite size behavior should be 
well reproduced by models in which both the external drift 
and the dissipation approach zero [5]. For these reasons, 
a new class of conservative sandpile models have been in- 
troduced, in which dissipation is prevented by periodic 
boundary conditions and no external input of grains is 
allowed. These models are called Fixed-Energy Sandpiles 



(FES) 6., since the total number of grains, or energy, is 
a constant of motion, fixed by the initial conditions. As 
a consequence, the activity of the system only depends 
on the total energy. For stochastic (i.e. Manna-like) up- 
dating rules, a threshold energy exists, above which the 
system does not relax to a stable state and a non-zero 
activity is dynamically maintained. Using finite-size scal- 
ing techniques, the critical behavior of the stochastic FES 
has been shown to belong to a particular class of absorb- 
ing state phase transitions (APT) , different to the directed 
percolation |7I8I9I1()| . On the contrary, statistical mechan- 
ics approaches miss to pinpoint the behavior of determin- 
istic FESs, with BTW updating rule gj. In fact, BTW de- 
terministic FESs (DFES) present very strong non-ergodic 
features due to the existence of periodic orbits in which 
the system eventually enters in all the range of possible 
values for the total energy |7I11| . In particular, in 
DFES has been studied on a square lattice with periodic 
boundary conditions (discrete torus) and it has been es- 
tablished that by varying the order parameter (the energy 
density) the system undergoes a transition from the frozen 
phase with an absorbing state to the active phase of even- 
tually periodic dynamics. Moreover, in the active phase, 
a "devil staircase" plateaus structure (strictly related to 
the behavior of the average periods) has been observed. 
The recently proposed exact solution for the dynamics of 
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the one-dimensional DFES ^2] corroborates the idea that 
deterministic BTW model couldn't be comprised into the 
same universality classes of stochastic or dissipative mod- 
els. 

The present work deals with DFES defined on an un- 
oriented connected graph, ft is focused on the relevance 
of dynamical invariants that are responsible for its non- 
ergodic features and undermine a purely statistical ap- 
proach. We note that a set <2>i, <P m of distinct constants 
of motion, ranging into vi,..,v m possible values respec- 
tively, determines a partition of the configuration space 
C into M = Y[" 1 Vk dynamically invariant classes of con- 
figurations (or "atoms" |23 m the language of dynamical 
systems). Such atoms are, in other terms, the counterim- 
ages of possible values assigned to the invariants. A natu- 
ral problem is then the determination of the maximal (i.e. 
the most refined) partition induced by such invariants. 

Dynamical invariants will be studied by recovering some 
important results already shown to hold for the dissipa- 
tive model In that case, with open boundaries and 
random addition of sand grains, invariants cannot play 
the same partitioning role. However, some of our prob- 
lems are implicitly shadowed there, and most of the alge- 
braic tools may be resumed and adapted to our model. In 
particular Ref. |13| suggests that dynamical toppling in- 
variants can be generated from the harmonic functions of 
the discrete Laplacian operator. The Laplacian can be de- 
fined on graphs by standard techniques of algebraic graph 
theory |14j . By means of modular algebra, we introduce 
a computational technique allowing for the evaluation of 
the harmonic functions even for quite large sizes. 

The group theoretical approach proposed by Dhar J3| 
has never been attempted in the conservative model. In- 
deed the Abelian Sandpile Group (ASG) is based on two 
dynamical properties: the addition operation and the exis- 
tence of a unique stable state for the subsequent relaxation 
(see Refs. |15llfi| for excellent reviews). In the conserva- 
tive model, there is no addition of grains and the final 
state, being a periodic orbit, is not univocally defined. A 
crucial point of the paper is the definition of the Toppling 
Group of "isoinvariant transformations" , that can be view 
as the natural extension of ASG for conservative models. 
Such an algebraic formulation, together with a relevant 
result of graph theory proves that the total number 
M{Q) of invariant atoms generated by modular harmonic 
functions is related to a remarkable topological property 
of the graph. Indeed N(G) equals the "graph complexity" . 
This quantity, defined as the number of the spanning trees 
of the graph, has been used in several contexts, such as 
the study of electrical networks going back to the semi- 
nal works by Kirchoff 18 .. The graph complexity can be 
evaluated by the product of the non-zero eigenvalues of the 
Laplacian matrix. More in general, the relation between 
physical properties and topology on graphs is a very in- 
teresting issue |19| . In the dissipative models, for instance, 
the use of conformal field theory has pointed out the rel- 
evance of boundary conditions in the computation of the 
critical exponents of correlation functions |2()j . 



Afterwords, for the relevant case of a L x L torus, [TT] 
we shall compute the complete class of independent in- 
variants for size L < 24 pointing out that the distribution 
of the invariants vs. L has a highly complex and unpre- 
dictable structure. In the meantime, the total number of 
invariant atoms has a very regular behavior which can be 
also analytically evaluated within the general framework 
of graph theory. In particular, the number of atoms grows 
exponentially with the size of the system. The large num- 
ber of toppling invariants provides a qualitative descrip- 
tion of the non ergodic behavior of the DFES. However 
such an explanation is not complete, since the basins of 
attraction (sets of configurations evolving into a periodic 
orbit or an absorbing state) are contained in the atoms, 
giving rise to a subpartition. By means of "isoinvariant 
transformations" we estimate the number of basins per 
atoms. Interestingly, such estimates provide a link with 
the "plateaus scenario" presented in JT]. The centers of 
the plateaus correspond indeed to peaks of the number of 
basins per atom. At least for small sizes, the calculated 
explicit form of the invariants allows for a direct explo- 
ration of the internal structure of the partition. Impos- 
ing the usual physical constraint of energy positivity at 
each site, we get that at low energies most of the atoms 
are empty; while at high energies all of them are filled 
by more or less the same number of configurations. As a 
consequence, a transition point E\ between these differ- 
ent regimes can be evaluated. We show that there exists a 
dynamical counterpart of such a transition, consisting in 
a qualitative change at E\ of the only relevant dynamical 
observables for a frozen system, i.e. the halting time and 
its fluctuations. 

The plan of the paper is the following. In Section [3 
we define the DFES model with BTW dynamics for a 
general unoriented graph Q. In Section |31 the class of in- 
dependent dynamical invariants generated by the discrete 
Laplacian is defined. Section rj] is devoted to the exten- 
sion to DFES of the algebraic approach of p^l by means 
of "isoinvariant transformation" and Toppling Group. In 
Sectional we compute the independent invariants for tori 
of size L < 24. A numerical study of the refinement of 
invariant atoms into basins is presented in Section [H] In 
Section^ a detailed study of the structure of atoms parti- 
tions is provided for small size (L = 4). Finally in Section^ 
the low energy transition induced by energy constraints is 
discussed. Conclusions and outlook for the future work are 
presented in Sectional The proofs of the theorems, a brief 
survey of modular algebra and the numerical algorithm 
used to find independent invariants are provided in the 
Appendices. 

2 The model: definition and properties 

Let us consider a generic connected unoriented graph Q 
with N vertices, labeled by integer i — 1, 2, . . . , N, con- 
nected pairwise by a set of unoriented edges defining a 
neighbouring relation i ~ j. At each site i, we introduce 
an integer variable z(i) £ Z, the number of sand grains in 
the site, physically interpreted as a local amount of energy. 
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A configuration Z is the observable Z : Q — > C, where the 
configuration space C — Z N is the iV-dimensional domain 
of integers. Moreover, each site is endowed with a fixed 
integer parameter, the critical energy z c (i) 7 such that if 
z(i) > z c (i), the site i becomes metacritical and di units 
of energy are equally redistributed among its neighboring 
sites, where the degree di of the vertex i is the number 
of its neighbours. Such redistribution event is called top- 
pling. The BTW dynamics consists of a parallel updating, 
in discrete time, of all metacritical states; therefore the 
evolution rule for Z t = {z t (i)} can be summarized as fol- 
lows, 

zt+ifl = zt® - Y^V%0[ztti) - z c (j)} , (1) 

3 

where 6{x) = 1(0) if x > (x < 0) and the integer N x N 
matrix V 2 is the Laplacian matrix |14| 

{di if j = i 
-1 if i~j (2) 
otherwise , 

In the following, the evolution rule (JIJ will be indicated 
with the operator Ubtw, i-e. Zt+i = UsTwZt- If the 
total energy E — J^i Z M i s sufficiently high, topplings 
propagate, creating avalanches of energy that cover the 
whole system. It is a common convention to assume the 
site variables bounded between and a maximal value 
imax = 2di — 1. The lower bound is simply due to the 
fact that negative amounts of grains are unphysical, while 
the upper one tends to avoid the useless treatment of en- 
ergy ranges forcing all sites to be active. Here, we consider 
the whole space Z N for a general approach, pointing out 
where necessary the effects of such additional constraints. 
The forward motion is completely deterministic, and the 
system eventually falls into a closed orbit, i.e. a fixed point 
or a limit cycle. Orbits {Zq, Z\, Z<i, ...} in C arc thus com- 
pletely determined by the initial condition Zq. 

An important remark concerns the symmetry imparl- 
ances of the dynamics. First the BTW toppling redistri- 
bution rule is completely symmetric, thus, given two con- 
figurations Z and Z' , if there exists a transformation £ 
(belonging to the symmetry group of the graph) such that 
£Z = Z' and £,~ 1 Z' = Z, then £ commutes with the dy- 
namics (£ oUbtw = Ubtw ° £)■ The second invariance 
is due to an internal symmetry of the evolution rule with 
respect to the transformation z(i) — ► z'(i) = 2di — 1 — z(i) 
Vi. 

Recently, Bagnoli and coworkers have investigated 
numerically the activity patterns of the DFES model on 
two-dimensional tori of different linear sizes L. Varying 
the energy density E = E/L 2 , the system undergoes a 
transition from a frozen phase (absorbing state) to the 
active phase characterized by eventually periodic dynam- 
ics. Moreover the phase diagram (the density of active , 
i.e. metacritical, sites vs. E), displays a plateaus structure 
organized as a devil-staircase In correspondence to 
these plateaus, the average length of the periods is small 
and approximately constant for increasing sizes L. A sim- 
ilar step-like structure with very short periods is partially 



recovered in the one-dimensional model, for which an ex- 
act solution is presented in Ref.|12|. Some of the intrigu- 
ing dynamical features of the two dimensional case will be 
qualitatively explained in Sectional 



3 Toppling invariants and discrete harmonicity 

The general concept of "toppling invariant" for a sandpile 
model on any kind of topology is quite obvious. Let Z t 
and Zt+i denote consecutive configurations along the or- 
bit started from Zq. A functional <P(Z) may be defined a 
toppling invariant if, for every t, 

${Z t+x ) =$(Z t ) . (3) 

The simplest conserved scalar form on Z N is the to- 
tal energy, that suggests to use linear integer functions 
of the site energies. On a periodic one-dimensional lattice 
the exact solution shows that constant of motions can be 
defined by means of linear scalar functions <P that are in- 
variant modulo the lattice size N . Such a property is a 
consequence of the simple linear periodic geometry; then, 
on a general graph, we should rather consider functions 
through a (mod-R")-linear form. 

On a generic connected unoriented graph Q (without 
multi-links and self-links), we introduce a set T{Q,K) of 
functions /:<?—> Zk defined on the nodes of the graph Q 
and taking values in the ring Zk = {0, 1, . . . , K—2, K — 1}. 
A KQ-functional $ } {K,Q,Z) generated by / G T(Q, K) 
is the quantity 

$ f (K, G,Z)=(j2 / W Z W j modif . (4) 

The range of a i^t/-functional is between and K — 1. 
Furthermore, the restriction f(i) S Zk is not relevant. 
Indeed, a substitution in Equation Qof / £ ^F(G, K) with 
/' : Q — > Z and f(i) — f'(i) mod K, gives rise exactly to 
the same functional. 

K 

Hereafter, we shall use the notation = for the equal- 
ities between elements of Zk- In such expressions, any 
integer has to be replaced with its if-modulo and all op- 
erations are to be intended in this "modular" sense. The 
space T(Q, K), endowed with the mod K sum, is a finite 
Abelian group, and therefore we will use the usual defi- 
nitions of finite group theory, such as element periodicity 
|22j . generator set, and group morphism. Relevant infor- 
mation on the i^-(Abelian)moduk) structure of J-(G, K) 
are summarized in Appendix 1X1 

Let us introduce the group morphism Vff : F{G, K) — » 
J-(G, K) which is the natural realization in J 7 (G, K) of 
the usual Laplacian operator defined by the matrix V 2 . 
In particular, for functions / G !F(G,K), V|-, is defined 

by 

(Vi/)(*):-(V 2 /)(z). (5) 
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A function h G J~(G, K) is said to be (discrete) K- 
harmonic on Q if 

(\7 2 h)(i)=0 VieG- (6) 

Alternatively, h is i^-harmonic on Q if it belongs to the 
kernel of Such Ker(V|-) is a subgroup of T{Q,K) 
and it will be denoted by TC(G, K). 

Functionals generated by -ftT-harmonic functions play 
a relevant role in FES dynamics. In particular, the KQ- 
functional ^/(K, Q, Z t ) is a toppling invariant if and only 
if its generating function / is i-T-harmonic in Q . The proof 
of this result is provided in Theorem 1 of Appendix [51 

The energy, generated by the constant harmonic func- 
tion c = 1, is the only i^-harmonic function for any value 
of K. Fixing the energy is equivalent to define the in- 
variant "energy surface" where the motion takes place. In 
analogy with the conserved quantities (motion invariants) 
of classical mechanics, Theorem 1 implies that the configu- 
ration space is divided into dynamically separated subsets, 
characterized by the value of ^/(K, Q, Z). In other terms, 
these invariants establish a partition of the energy surface, 
i.e. an exhaustive collection of disjoint subsets, that in the 
dynamical systems theory are called atoms It is wor- 
thy noting that not all the i-T-harmonic functions are inde- 
pendent. For example, for f,g£ H(G, K), if / and g = 2f 

then <P g (K, G, Z) = 2<P f (K, G, Z). Therefore, for each con- 
figuration Z, the value of <P g (K,G , Z) is determined by 
<P f (K,G,Z). 

In general, the function h G H(G, K) is defined to be 
dependent on the functions h n G 7i.(G, K n ) if there exists 
a set of numbers aj such that for any function Z, one has 

N 

$ h (K, G, Z) = a E(Z) + a^ hj (Kj,G, Z) (7) 

3=1 

The energy of the system E(Z) has been considered apart, 
so that the independent invariants can be defined neglect- 
ing an arbitrary constant. 

For instance, in a one-dimensional system the enu- 
meration of all the independent functions is trivial; the 
structure of the discrete Laplacian imposes N — 2 con- 
straints on N variables. The two independent functions are 
for instance those corresponding to the total energy and 
the "linear momentum" mod-/V along the periodic system 
(h,2(i) = z, Vz) . Increasing the dimensionality of the lattice 
or changing the topology to a generic graph, the structure 
of the configuration space becomes more and more com- 
plicated and the number of dynamical invariants rapidly 
grows with the system size. 

In general, a complete set of independent dynamical 
invariants can be computed by the Smith decomposition 
of the Laplacian matrix V 2 (see Ref. ^SJ) as 

S = A\7 2 B = diag(0,g 1 ,g 2 ,... 1 g N _ 1 ) , (8) 

where A and B are integer matrices with determinant ±1 
(unimodolar matrices) The k— th column of B in JHJ de- 
fines an independent harmonic function mod gk ■ We have 
two important corollaries of this result: 



1 - the total number of independent invariants is finite; 

2 - the number of atoms defined by such non constant 

invariants is M{G) = rifcTi 1 9k- 

An explicit calculation via Smith decomposition is pos- 
sible only for graphs of small sizes, because, for intrinsic 
features of the algorithm, integer numbers greater than 
the maximum allowed by the computer are early involved. 
For example, on L x L tori, the transformation JSJ) may 
be worked out only up to L = 5. An alternative approach, 
described in details in Appendices C and D, allows to com- 
pute the invariants for tori with L up to 24. Unfortunately, 
this alternative method cannot assure the completeness of 
the set of invariants. In our experiments, completeness is 
guaranteed only for K < Kg = 20000. This is once again a 
computational bound. More in general, for a graph G and 
a suitable integer Kg the algorithm provides a set Ig.K g 
of independent functions hj g HiG^j 3 ) of periodicity 
(pj are primes, Lj are integers and < Kg Vj). Further- 
more, every function h € H(G,K) (K — Yij=T ®j w ith 
0j E N, 6j < Kg) depends on the set Ig,K g - 



4 Toppling Group and Isoinvariant 
Transformations 

On the configuration space it is possible to introduce a 
class of transformations preserving the values assumed by 
the constants of motion. In particular, an "isoinvariant 
transformation" is a mapping r\: Z G C — > Z' G C such 
that <P h (K,G,Z) = <P h (K,G,Z') for any integer K and 
any if-harmonic function h. 

Theorem 2, reported in Appendix iBl states that a nec- 
essary and sufficient condition for a transformation to be 
isoinvariant is that it can be written in the form 

Z' = rj(Z) = Z + V 2 U (9) 

where U G Z N is an integer vector. 

We observe that the class {ri} iso of transformations, 
defined by 0, operating linearly on the configuration 
space, constitutes an Abelian group. A set of generators 
for this group are the elementary toppling operators 7]i, 
decreasing by di grains the site i, and augmenting all its 
neighbors by 1 unit. Hence such group will be called the 
Toppling Group. The BTW dynamics itself is an example 
of isoinvariant transformation, depending however at each 
step on the metacritical domain of the state. This makes 
in general t/UbtwZ UbtwV^- (Of course, in a reduced 
configuration space C, we should admit only those trans- 
formations ensuring that Z' belongs to it). 

Some global properties of the toppling group can be 
studied by means of algebraic graph theory (see for in- 
stance Refs. Particularly important is the concept 
of preflow that has been introduced in |17j . 

Given a graph G, we call integer "preflow" an integer- 
valued function defined on its edges after having assigned 
them an orientation. The preflows constitutes a group with 
the addition in Z. 
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The isoinvariant transformations are a subgroup of the 
preflows on Q. In our context, we can apply a relevant re- 
sult on integer preflows stated in Ref. [T2|. More precisely, 
in Appendix [E] we will prove that the number Af(G) of 
atoms generated by toppling invariants of the form 
equals the complexity of the graph Q, i.e. the number of 
its spanning trees [14| . For example, in a generic tree- 
graph the number of spanning trees is 1 , while on a chain 
is the number of sites. Moreover, this number can be eval- 
uated from the Laplacian matrix V 2 as the product of all 
non-zero eigenvalues , i.e. 

1 N 

W) = jvrII A *- ( 10 ) 

i=2 

(It is a classical result of algebraic graph theory that the 
quantity (fTTTfi is an integer). 

A sort of partition into atoms and isoinvariant trans- 
formations, however with different meaning and physical 
interpretations, were already defined in Dhar's algebraic 
approach |13I15| . In the dissipative model, partitions do 
not correspond to constants of motion (the system remains 
within an atom during an avalanche, jumping to different 
atoms by addition of new grains). Another relevant prop- 
erty of dissipative model is that in each atom there is one 
and only one recurrent configuration, this property does 
not have a counterpart in FES, where recurrent config- 
urations cannot be defined. Moreover, from an algebraic 
point of view, we note that in the open case, in order to 
take into account dissipation, the matrix playing the role 
of the Laplacian V 2 is not singular, permitting a set of 
operations not allowed in our case (e.g. the definition of 
atoms using the inverse of the matrix). Anyway, in the 
link with the graph complexity is also put into evidence 
following a different path, the so called "burning test", 
instead of our preflows criterion. 
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Table 1. Number Q p ^ of independent harmonic functions of 
periodicity p L for L between 2 and 13 and Ki = 20000. 



that no simple rule can be inferred to foresight results 
corresponding to larger sizes. On the contrary, the total 
number of atoms M{G) seems to have a quite regular be- 
havior. Indeed, Af(G) can be evaluated as 

wo) = n p" QpL ■ (n) 

p L <Kg 



5 Invariants on a torus 

In this Section we focus on the case in which Q is a two- 
dimensional regular periodic lattice of size N = L 2 , i.e. a 
torus; in this system, the non-ergodic properties of FES 
dynamics was observed for the first time in jolllj . The 
toppling invariants partition will be studied by means of 
the algorithm introduced in the appendices IU1 and iDl By 
considering Kg — 20000 (the bound Kg is introduced in 
Section[2Jl, we are able to compute the invariants for sizes 
L up to 24. The results for L < 13 are reported in Tabled 
For such sizes the algorithm seems to be efficient to find a 
complete set of invariants; indeed, the largest periodicity 
appearing in Table HI is much smaller than Kg = 20000. 
For L < 5 this has been confirmed by a direct evaluation 
of the whole set of independent invariants by means of a 
Smith decomposition (as already pointed out, the Smith 
decomposition cannot be easily applied to larger systems). 
In Tabled we call Q p l the number of functions of period- 
icity p L in the class Ig t K g of invariants. 

Table evidences that the number of invariant as 
functions of L, p L and Q p ^ is very irregular and it seems 



We exploit the property that an invariant generated by 
a .ftf-periodic function can have K values, and we assume 
that Ig,K g is a complete set. The logarithmic plot of N(G) 
vs. L reveals a clearly quadratic behavior (see Figure [SJ). 
More precisely 



Af(G) ~ exp(cL 2 ) 



(12) 



where c = 1.20 ± 0.05. In the figure all L < 24 and some 
of the larger sizes have been considered. For some values 
of L we get a different behavior, however it is likely that 
such anomalies could depend on invariants of periodicity 
larger than Kg = 20000, which are not captured by the 
actual implementation of the algorithm. M(G) can also be 
evaluated by means of formula l|10f) in this case the non- 
zero eigenvalues of the Laplacian matrix are analytically 
known, this gives 



L-l L-l 



^)=^nn 



1=1 m=l 



4 - 2 cos 



2tt£ 



2 cos 



27rm 



(13) 
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Fig. 1. The natural logarithm of the total number of atoms 
log(7V(<7)) as a function of the torus size L. 



For large L, AT(G) ~ exp(ct^L 2 ), where the parameter c t h 
can be computed as [T5] 



Cth 



dad/3 
~4^~ 



log (4 - 2 cos a- 2 cos 0) ~ 1.17, 

(14) 



that confirms the results obtained by means of numerical 
evaluations. 



6 Atoms and basins 



Number of 
basins 


Orbit period I 


Basins of period I 
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16 
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418 
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Table 2. Some typical atoms partitioning in different basins 
of attractions at the energy density E — 2.25. The first column 
reports the number of basins contained in the atom. The next 
columns display how such number is distributed for basins cor- 
responding to orbits of different period length. The data are 
obtained numerically averaging over 200 atoms and 1000 initial 
conditions for each atom. 



The exponential growth of the atoms for the DFES on the 
torus provides a qualitative explanation of its non ergodic 
features, since it implies a fragmentation of the config- 
uration space into a very large number of dynamically 
intransitive domains. However, such an invariants-based 
partition does not resolve completely the non ergodicity. 
Dynamical basins of attraction determine indeed a further 
non trivial subpartition, whose properties will be investi- 
gated in this section by means of isoinvariant transfor- 
mations introduced in Sec. 4. By iterating such transfor- 
mations, one can generate indeed a sufficient number of 
states for a reliable estimate of the basin weights within 
every single atom, revealing its internal structure. 

In particular, for a torus of size L = 40, at differ- 
ent energy densities, we collected data for 200 atoms and, 
in each atom, 1000 initial states (obtained by extracting 
random integer vectors for the transformation ©)■ Some 
typical results for an energy density E = 2.25 are dis- 
played in Table J5J. Each couple of rows corresponds to 
a different atom. The first column reports the number of 
different basins identified over 1000 initial configurations. 
In the next columns basins are grouped according to the 
period of the final orbit (in the first row there is the orbit 
length, in the second the number of configurations evolv- 
ing into an orbit of such length). The energy E = 2.25 
corresponds to a plateau where the average orbit length 
is eight (see ^D). Nevertheless, due to finite size effects, 
a wide range of orbit lengths is present. Different behav- 
iors can be observed. In particular: 1) there are atoms 
divided into very few (or even one) basins; 2) there exist 




E 



Fig. 2. The average number of different basins, emerging from 
1000 states of the same atom, vs. E. In the average, 200 atoms 
are considered for each energy. This behavior does not depend 
on the number of starting configurations or the lattice size. 



atoms divided into many (up to hundreds) basins. A rea- 
sonable property (frequently observed indeed) is that the 
former case corresponds to very long orbits, and the latter 
to short ones. However, there are remarkable exceptions. 
It is also noteworthy that, when an atom is divided into 
many basins, only few periods are present, and all of them 
are simple multiple of a minimal period. This is true for 
all energies. 
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In Figure |3 we plot, as a function of E, the average 
number of distinct basins reached starting from 1000 con- 
figurations belonging to the same atom (the average is 
obtained considering 200 different atoms at the same en- 
ergy). We checked that this behavior does not depend on 
the number of starting configurations or the lattice size. 
Indeed, by varying such parameters, we obtain an analo- 
gous figure after a suitable rescaling of E and B. A com- 
parison with data in Ref. shows that the plateaus char- 
acterized by short periods correspond to peaks in Figurc|3 
Hence, in the middle of the plateaus the atoms are highly 
fragmented into many basins, while in the transition re- 
gions between different plateaus there are few basins per 
atom. This suggests that the dynamics of the system plays 
a non trivial role in the refinement of invariant atoms, con- 
firming that there is a structural difference with respect 
to the dissipative case: indeed, in the Dhar's algebraic ap- 
proach, the partition into atoms provides a complete de- 
scription of the configuration space, since to each atom 
corresponds a single absorbing (recurrent) state (see Sec- 
tio |2J . The possible existence of further simmctrics in- 
duced by the deterministic BTW rule suggests that the 
DFES model belongs to a different class of universality. 
Similar simmctry based arguments may be used to distin- 
guish the DFES model from stochastic models. 



7 Atoms filling for two-dimensional lattices 
of small size 

Up to now we assumed that the configuration space C is 
extended to the whole Z N . However, for reasons recalled 
in Section [21 in sandpile models the energy is generally 
assumed to be non negative, and it is also bounded by a 
maximum value. Hereafter, we assume that on a toroidal 
lattice < < 7 Vi, calling this reduced configuration 
space C C C. 

A natural problem arising from such constraints is the 
way atoms are filled by configurations. Indeed, while in C 
there are infinite configurations compatible with the in- 
variants, in C their number is finite, and it depends on 
E. Since this number exponentially grows with the size 
N = L x L, an exhaustive numerical analysis is bounded 
to L = 4 (see Table for a complete set of invariants in 
the matricial representation). 

Let p(rn) be the number of atoms containing m al- 
lowed configurations. The average occupation m and its 
standard deviation are then given by: 



Em mp(m) 

T, m p( m ) 



Am 



ErJrrc - m) 2 p{m) 
Em P(m) 



(15) 



The exhaustive calculation for L = 4 provides p(m) at 
different energies. The results are plotted in Figures and 
At very low energies, most of the atoms are empty. In 
this case they will be denoted as virtual atoms, since no 
permitted configuration exists realizing the corresponding 
values of the invariants. The upper panel of Figure|3|shows 
that for small energies p(m) is maximum at m = 0, then 
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Table 3. Complete set of independent generators for the non 
constant harmonic functions for a two-dimensional FES of size 
L = 4. 



E=10 

* E=12 

o E=14 

+ E=16 



Fig. 3. The logarithm p(m) versus m (p(m). The upper panel 
is relevant to low energies E where p(0) is maximum and an ex- 
ponential decrease is present. The lower panel refers to higher 
energies where the occupation is peaked around the average 
value fa. 



it exponentially decreases with m. On the other hand, 
the lower panel evidences that for higher energies p{m) 
vanishes at m = 0, so that all atoms are filled by allowed 
configurations. In particular, p(m) is peaked around the 
average occupation m. Finally, we note that rh rapidly 
increases with the energy of the systems. 

Let us now analyze the shapes of the peaks. In the up- 
per panel of Figure0J data calculated for E — 28 are fitted 
using a Gaussian curve with parameters corresponding to 
the actual average value m and standard deviation Am. 
p(m) has significant deviations from a Gaussian behavior, 
with an asymmetric shape exhibiting a long tail for large 
occupations. The lower panel displays the average occupa- 
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Fig. 4. In the upper panel a detailed plot of p(m) for E = 28. 
The comparison with the Gaussian distribution (same average 
and standard deviation) shows an asymmetric behavior. In the 
lower panel the standard deviation Am vs. the average occu- 
pation number m. The comparison with the curve m 1//2 shows 
that atoms are not filled randomly. 



tion m versus the corresponding standard deviations Am 
for different energies. A random occupation would imply 
Am — \fm~, while the figure points out a deviation from 
such behavior. Precisely, the standard deviations are much 
larger than y/Wi, a signal of correlation in the occupation 
number. 

Since we imposed the constraint < z% < 7, an anal- 
ogous transition is present also for E ~ 7 because of the 
symmetry mentioned in Section [21 It would be important 
to establish the robustness of this scenario at higher L, but 
direct numerical experiments are computationally very de- 
manding. However, for some of the properties, we have 
checked that our results are consistent also at larger sizes. 

8 Low energy transition 

Assuming that previous results hold independently of the 
system's size and topology, there should exist a transition 
induced by the deformation of the space C with respect to 
C = 7i N , affecting only low- (high) energy hyper-surfaces. 
To identify the energy of such a crossover, we need to 
evaluate C(E, N), i.e. the number of configurations having 
energy density £ on a system of size N. Let rii be the 
number of sites having energy i (i = 0, 1, ... , i m ax) and 
Xi = rii/N; then the number of configurations C(E, N) is 

C(E,N) = J ^L-dxi~ J e- N T Ji ^^) dXi (16) 

In a saddle-point approach, for large N, the values of 
are evaluated by minimizing ^ Xi lqg(a;j), with the con- 
straints . Xi = 1 and Yli = By means of La- 
grangian multipliers, the minimum of Xi log(xj) can 
be easily obtained in terms of E and i ram . 

C(E,L) ~ e b< — . (17) 

Where bi max (E) is a function that can be analytically eval- 
uated. By comparing (| 1 3|1 and <|17(1 . whenever bi ma:c (E) < 
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Fig. 5. The average halting time (solid line) and its fluctua- 
tions (dashed line) vs. the energy density (fluctuation are mul- 
tiplied by 10, in order to put plots in the same figure). The 
vertical dotted lines denote the crossover density energy Ei. 



Cth the number of configurations results to be exponen- 
tially smaller than the number of atoms, and most of them 
do not contain any configuration. On the other hand, for 
bi m a X (E) > c th, all possible atoms are expected to be occu- 
pied. Therefore, the crossover energy E\ can be obtained 
as bi max (E) = Cth, which gives E\ ~ 0.72. 

Now we want to check if such a crossover has also dy- 
namical consequences. Since at E\ the system is frozen 
(indeed periodic orbits are present only for E > E ~ 2.06 
jll|). the only relevant dynamical quantities are the av- 
erage halting (freezing) time r and its fluctuations At. 
Now, let tz be the halting time for the orbit starting from 
a configuration Z. We define 

= Ez(E) T Z(E) = Ez(E) T Z(E) _ t 2 flg x 

N E N e ' y J 

where the sums run over the configurations Z at energy 
E, Ne being the number of such configurations. In the nu- 
merical simulations we have averaged over 5000 different 
initial conditions. 

In Figure t and At are plotted versus the energy 
density. Since, for such low energies, the halting time is 
small, we can perform simulations for very large systems 
(here, L — 800). The vertical dotted line marks the tran- 
sition between the regions where the number of atoms is 
smaller or larger than the number of configurations. In the 
former region, r is characterized by a series of slow-downs 
forming small plateaus where r results to be an integer 
value. The minima of fluctuations occur in correspondence 
of these plateaus, while their maxima are in the transitions 
between different plateaus. Such behaviors of r and At 
can be directly explained if for an energy corresponding 
to a plateau almost all of the configurations are character- 
ized by the same halting time. For energies between two 
different plateaus, only the two halting times characteriz- 
ing the near plateaus are possible. If the dynamics freezes 
with equal probability in n or in n+1 steps (all other halt- 
ing times being not allowed), then At = 1/4, indeed the 
height of the fluctuation peaks results to be about 0.25. 
On the other hand, at energies larger than E\, the be- 
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Fig. 6. The fluctuations of the halting time for different lattice 
sizes. The vertical dotted lines signal E\. Below E\, fluctua- 
tions are oscillating functions of the energy density. For E > E\ 
fluctuations increase with the energy and seem to be indepen- 
dent of the lattice size, corroborating the conjecture about the 
existence of a transition at E\. 



haviors of r and At are much more regular. This should 
follow from the large occupations of the atoms, permitting 
many halting times at the same energy. In such case, the 
system cannot be described by the simple picture used for 
E<E X . 

In figure El we plot the fluctuations of the halting time 
as a function of the energy density for different torus sizes. 
In the low energy region, by increasing the size of the 
system oscillations become more and more rapid and the 
peaks move towards E = 0. This is consistent with our 
picture, since at a given low energy density, the probability 
that the evolution freezes into n steps increases with the 
size. Indeed, this probability depends on the number of 
way a configuration can contain one or more blocks of sites 
that freeze in a certain time. On the contrary, above the 
critical energy E\ , the fluctuations seem to be independent 
of L, which is a rather surprising result. 



9 Conclusions 

In the present work we have investigated the structure of 
the configuration space of DFES models on generic unori- 
ented graphs. Our major result, the complete and algorith- 
mically explicit calculation of the toppling invariants, ex- 
tends to the conservative case the group theoretical frame- 
work introduced by Dhar for dissipative sandpilcs. The 
Toppling Group seems to be a very general feature link- 
ing sandpiles dynamics with graph's algebraic properties. 
Using this algebraic approach it is possible to identify the 
exact way the configuration space is partitioned by the 
dynamics into invariant subsets, and to determine their 
main properties. 

The validity of the analytical results is corroborated 
by an independent numerical analysis carried out in the 
case of two-dimensional lattices. 

As by-products of our analysis, many properties of the 
two-dimensional system are elucidated. In particular we 



give a qualitative explanation for the abundance of or- 
bits with very short periods. In addition, the further re- 
finement of atoms into attraction's basins reveals an un- 
expected relation with the devil's staircase structure of 
Ref. . The centres of the plateaus correspond to en- 
ergies displaying peaks in the number of basins per atom. 
We argue that the absence of a one-to-one correspondence 
between periodic orbits and invariant subsets makes the 
deterministic conservative model "structurally" different 
from both dissipative ones and the stochastic conserva- 
tive models, supporting the thesis of a distinct universality 
class for the deterministic FES. 

We also point out that any further constraint can have 
relevant consequences on the organization of the configu- 
ration space into invariant sets. For example, the usual 
constraint on the energy positiveness introduces a transi- 
tion in the frozen phase, saparating the region where most 
of the atoms are empty from the region where atoms are 
filled with a growing number of configurations. For small 
size systems, this analysis can be performed by an ex- 
haustive counting procedure. On the other hand, for large 
systems, the transition energy may be analitically esti- 
mated by asymptotic methods. This transition has also a 
dynamical counterpart in the behavior of the halting time. 

This work may as well be the basis for a more general 
study on the the statistical properties (spectral features, 
noise etc.) of stochastic counterparts as dynamical systems 
presenting a mixture of randomness and regularity due to 
the existence of a partition of the configuration space with 
an external weak random perturbation (see |24j for other 
examples) . 
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A Modular algebra 

Since the definition of i^tj-functionals involves the mod- 
ulo operation, we introduce some useful definitions and 
notations of modular algebra (for details see e.g. 

Definition 1 The functions ofT{Q,K) are naturally en- 
dowed with a structure of K-(Abelian) module (i.e. a vec- 
tor space where scalars belong to the ring *Lk and not to 
a field). In particular, the sum of f and g is denoted as 

h= f + g and the product of f and a 6 Xk as h = af . 
is the neutral element of T(Q,K). 

Let A4(K) be a generic K- module. An element m £ 

M(K) has periodicity P e Z K , P < K, if Pm = 0, in this 
case P is a factor of K. We shall now proceed in extending 
as far as possible the basic notion of linear algebra in the 
modular sense. 

Definition 2 A function F : M(K) — > M{K) is a K- 
module homomorphism if for any m, n S M.{K ) and a,b € 

Zjf F(am + bn) = aF(m) + bF(n) . The kernel of F is the 
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set Ker(F) = {m £ M(K)\F(m)=Q}. Ker(F) is a sub- 
module of A4(K). The discrete Laplacian V^- in Eq. is 
an example of K -module homomorphism V^- : T{L, K) — > 
F{L,K). 

Definition 3 Consider a set S — {s a } C A4(K), with 
a = 1, . . . , N. The elements {s a } are algebraically depen- 
dent if there exists s@ £ S and a set of a a £ such 
that 

si3 = ^a a s a (19) 

Otherwise the elements of S are algebraically independent. 
Let Q C M. (K) be a set of algebraically independent ele- 
ments {g a }- The element ofQ are independent generators, 
if for any m £ M. [K) there exists a set a a £ Zk such that: 

m = ^2a a g a (20) 

ot 

A subset E C Ai(K) of independent generators {e a } is a 
basis if 

Y,a a e a = (21) 

a 

implies a a = for all a. If M(K) has a basis than M(K) 
is called a free K -module. We point out that any K -module 
has a set of independent generators, while there exist mod- 
ules where it is not possible to introduce a basis. 

Free iT-(Abelian)modules have very peculiar proper- 
ties, and they are similar for many aspects to the usual 
vector spaces. In particular any element m £ Ai{K) can 
be written in a unique way as (a a £ Z#): 

m = ^2a a e a (22) 

n 

As a consequence, in a free iT-module any homomorphisms 
F can be represented by means of matrices. Let us expand 
the generic element m £ A4(K) in the basis {e a } as in 
Ij22(l . The components b\, 62, ... of F{m) in the basis {e a } 
result to be: 

bp = >^ F/3 ja a a (23) 

a 

where the matrix Fp. a as in the case of standard vector 
spaces is defined by 

F(e a ) = J2 F P,c*ep ( 24 ) 


J-(Q ', K) is a free module and the homomorphism F : 
T{L, K) — > J-(L, K) is naturally represented as a matrix 
in the basis e m (i) — Si, m . However, the existence of a basis 
for a generic A'-module is not guaranteed (the if -module 
is not free). In particular, some submodules of T(Q, K) are 
not free. On a torus, let us call x the submodule generated 
by: 

91= (Jo) ^=(20) (25) 



i.e. any element of c £ x is a combination of the form 

4 

c = a%gi + d2(?2- X is Zzi-module with g\ and gi as inde- 
pendent generators. However x does not admit any basis. 
Hence, a representation of the homomorphisms Fx — ► X 
in matrix form is not possible. We note that the element 
<?2 bas period P = 2 < K = 4. 

B Theorems 1 and 2 

In this section, we provide a simple proof of Theorems 1 
and 2 for a generic undirected graph Q with N vertices. 
Let us define the sets 

R k = {i : z(i) = k} , (26) 

we call critical and metacritical regions R c and R m re- 
spectively 

R c = {i : z(i) = z c (i)} , and R m = {i : z(i) > z c (i)} . 
The toppling matrix T t is defined by 

T t = Z t+1 - Z t . (28) 

Theorem 1 The KQ -functional <Pf(K, Q,Z t ) is a toppling 
invariant if and only if its generating function / is in- 
harmonic on Q. 

Proof First we show that if-harmonicity is a sufficient 
condition. Let R m {t) be the metacritical set R m at time 
t indented with its non metacritical neighbors. Let i?„(i) 
be the complementary set, i.e. R m (t) U Rm(t) — Q. At 
time t, the sum defining •Pf(Z t ) in J3J may be split in two 
sums running over R m {t) and R^ n (t) respectively: 

M=Y,tif(i) zt(i)i£Rrn{t) (29) 
Bt=Y^ =1 f^zt{^)^£R' m {t) (30) 
Similarly we define: 

At+i =E£i/(») *+i(0 * e R m (t) (31) 
Bt+i =E£i /(O *t+i(0 * e KM (32) 
Since i?„(i) , is not touched by evolution, B t +\ = B t and 

$ f (Z t+1 )=A t+1 +B t+1 =A t+1 + B t (33) 
At+i may be rewritten as 

N 

A t+1 ^f(i)[z t (i)+T t (i)], (34) 
i=i 

for i £ R m (t). Consider now the difference 

N 

A t =A t+1 -A t = J2f(i)T t (i), i£R m (t). (35) 

i=l 
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The Abelian property of the BTW-rule ensures that the 
sum defining A t , can be rewritten as a sum over the 
strictly metacritical set, by writing: 

4 t = £[-*/(») + £/Cj)] 

«'=1 \ 3~i ) 

N 

i=l 

where we use the hypothesis of /c-harmonicity of /. There- 
fore 

$ f (Z t+1 )=A t +B t + A t =A t +B t = $ f (Z t ) (37) 

leading to the conclusion that _ftT-harmonicity implies top- 
pling invariance. 

Now we prove that if-harmonicity is a necessary con- 
dition. The hypothesis is that, for every Z t 

$ f (Z t+1 )^<P f (Z t ) (38) 

and the claim is (V 2 /)(m) = 0, for an arbitrary m. Let in 
particular Z t (i) = cWi, m , then, Z t+ i(i) = J2j~ m Tnis 
means 

# / (^ t+1 )-%(Z t )^(V 2 /)(m) = ) (39) 

□. 

Theorem 2 A necessary and sufficient condition for a 
transformation to be isoinvariant is that it can be written 
in the form 

Z' = n(Z) = Z + V 2 U (40) 
where U G Z N is an integer vector. 

Proof Two configuration related by formula © are clearly 
isoinvariant since 

N 

$ h (K, g, z') £ Y, h ^ ( z « + ( v2[/ )(*)) 

N N 
i=l i=l 

= <P h (K,G,Z) (41) 

where we used the symmetry properties of V 2 and the fact 
that h is a if-harmonic function. 

Let W = Z' — Z be the difference between the configu- 
rations Z and Z' belonging to the same atom, we have that 
B^W — SL, where B and S are the unimodolar and the 
diagonal matrices appearing in the Smith decomposition 
iJHJ, and L is an integer vector. The Smith decomposition 
yields B^W = B^\J 2 A^L and then W = \7 2 A^L, therefore 
W can be obtained applying the Laplacian matrix to the 
integer vector A*L □. 



C Independence in H(Q,K) 

The algebraic independence (|19|) between functions of the 
if-module Tt(G, K) implies that they are also indepen- 
dent according to definition Q. This is a simple conse- 
quence of the fact that the toppling invariants are "lin- 
ear functional" in the space J-(G,K). Therefore the par- 
tition of the energy surface generated by the whole set of 
if-harmonic functions H(G, K) is the same as the parti- 
tion generated by a subset of independent generators of 
H(Q, K). On the contrary of !F(Q, K), H{G, K) is not in 
general a free module l|21l) . In particular there may be 
generators of periodicity P < K . 

Let us now consider the harmonic functions correspond- 
ing to different if's. Since only the constant functions be- 
long to TC(G, K) for any K, the energy is a special invari- 
ant to be considered apart, as in definition 10. To this 
purpose we denote with H'(G, K) any set of generators of 
H(G, K) such that c(i) = l,Vz belongs to H'(G,K). The 
following theorem shows that the partition of the energy 
surface generated by all the if -harmonic functions can be 
evaluated by considering only K's of the form K — p K , 
where p is a prime number and re is an integer depending 
on p, i.e. considering the invariant generated by the func- 
tions of H '(GtP* 1 )- Let us first introduce a lemma which is 
a basic properties of integer number: 

Lemma 1 Let K = p^p^ 2 ■ ■ ■ Pr r the decomposition of 
K into prime factors p n . Any element of b G Zk can be 
written in a unique way as: 

r 

&=E«» 6 « (42) 
n=l 

where b n G Zj p >««i = {0, 1, . . . ,p* n — 1}; and the integers 
q n are given by: 

1n= II Ph ( 43 ) 

Ph¥=Pn 

Theorem 3 Let us consider the module H.(G,K), and let 
K = p^p^ 2 ■ ■ -Pr r the decomposition of K into prime fac- 
tors p n - For any h € TL{G,K) there is a set of functions 
h n S 'H{G,Pn n ) suc h that h depends on the h n 's according 
to definition Q). 

Proof ^From Lemma 1, any elements h(i) of the function 
(matrix in two dimensions) h G Ti(G,K) can be decom- 
posed in a unique way as: 

r 

h(i) = J2qnh n (i) (44) 

71=1 

where h n (i) G Z{ pK „} = {0, 1, . . . ,pn - 1}. Equation (|4*4"jl 
defines the functions h^s. Indeed we have 

r 

0K = Vk(ft) = £>»Vi-(/i») (45) 

n=l 

Lemma 1 entails that the neutral element K of H(G, K) 
can be obtained, by means of factorization (|44(1 , only as a 
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combination of the neutral elements Q{ p ~>»} of TC(G,Pn n )- 
Therefore, Equation l|45|) means that, for all n, h n £ 
T~i(G,Pn n )- Finally for the if^-functionals we have 

r / n \ 

n=l Vi=l / 
r 

= E"A»K"-^ Z ) ( 46 ) 

71=1 

in the last expression we used the equality q ■ a = q ■ 
(a mod (K/q)) which holds when K/q is an integer, (here, 
from K/q = p% n ). □. 

We now introduce a theorem which allows to simplify 
the study of H(G,p') for different values of the integer i. 
The following Lemma can be directly proved resorting to 
the properties of modular spaces. 

Lemma 2 A periodic element h £ TL(Q,p K ) of periodicity 

p L , i.e. p L h(i) = 0, can be written as h(i) — p K ~ L h'(i) with 
h'£ H{G,p l ). 

Theorem 4 Let us consider an integer k such that the 
any function h £ Ti(G,p R ) can be written as h(i) — h'(i) + 
c where c is a constant and h' is a function of periodicity 
smaller than p K . Any function belonging to TL{G,p e ), for 
any integer e, is dependent on a function h £ Ti(Q,p K ). 

Proof A generic function belonging to TL(G,p e ) will be 
denoted as ftw. First we focus on the case e < k. 
For any element h^ e ' £ TL(G,p e ) we consider the ele- 
ment h' £ J 7 (G,p K ), given by h'(i) — p K ~ t h(i). Since 

(V 2 /i')(«) =p h - £ (V 2 /i)(i) =p h - e (a p £ )=0 (a is an inte- 
ger), we have h' £ H{G,p h ). Moreover, <P h(e) (p e ,G, Z) = 
P € ~ K @h'(p K , G, Z) and this concludes the first part of the 
proof. 

We pass now to the case e > k. If all the functions 
ftW £ U{G , p K ) are a combination of a constant and a 
function of periodicity smaller than p K , Lemma 2 implies 
that 

=ph {K - l) {i)+c (47) 

with h^ k ~^ £ H{G 1 p K ~ 1 ). First we show by induction 
that in any module Ti.(G,p e ) with e > k all the functions 
h^' are a combination of a constant function and of a 
periodic harmonic function of period smaller than p K , i.e. 
from Lemma 2 

(») = p'-^h^-V (i) + b (48) 

In the case e — k this property is verified by hy- 
pothesis. We now suppose that the property holds for 
M t_1 ) and show that this is true also for h L . We de- 
fine h*(i) =p l ~ K h( L \i). Since p K h* = 0, from Lemma 2 we 
have that h* = p L ~ K h^ with £ H(G,p K ). We have 

p L - K hM\i) = h*(i)=p L ^ K h ( - K ^(i) and multiplying both 



sides by p K ~ x we get p' -1 /iM(i) = p L 1 h ( - K \i). Now h (K) 
satisfies <|47ll and therefore 

p'- 1 h^{i) p =p L - 1 {ph {K -V(i) + c) p =p'- 1 c (49) 

Equation Q49JI proves that ftW(i) — c is periodic of pe- 
riod p^ 1 and then from Lemma 2 ftw(i) — c = p/)/ t_1 )(z) 
where h^-^{i) £ HiG^ 1 )- Since by induction we 
supposed that Equation (|48|l holds when e = t — 1, 
we get ftW(i) - c = p(p u - K h {k - 1 \i) + b) i.e. h^(i) = 
p L ~ K+1 h( k ~ 1 \i) + d where d is an integer constant. This 
concludes the prove by induction. Finally, we consider 
the toppling invariants <Ph>{p e ,G, Z), from l|48|) we have 

<P h { P \ g, z) p =p^+ 1 ^-, ) ( P K -\G, z) + dE{z). □ 

We note that Theorem 4 does not prove that the inte- 
ger k exists. This can be proved, for instance by means of 
Smith decomposition (Section 0. 



D Algorithm 

D.l Matrix equations 

Let us consider a free module Ai(K) and its basis e a . 
Since any homomorphism F : Ai(K) — > Ai(K) is repre- 
sented by a matrix Fp >a , the equations determining the 
components k a of the elements of Ker(F) are 

M 

J2 F P,<*k*=0 (50) 

a=l 

where M is the number of elements of the basis. Let 
us show that the system (JSUJ, can be faced by resorting 
to techniques similar to those usually adopted in ordi- 
nary vector spaces. First, the following operations do not 
change the number of solutions of l|50|l . 

— Substitution. If the clement -F/3, 7 of the matrix defining 
system (|5U|I equals to 1, the corresponding element k 7 
can be eliminated by the substitution 

k 7 = - ^ F P,<* ka ( 51 ) 

Equation 1)51(1 can be discarded reducing the system 
from M to M — 1 equations and variables. 

— Columns and rows exchange. It is possible to swap two 
rows or two columns of the matrix Fp }a . This corre- 
sponds to a relabeling of equations and variables. 

— Row sum and difference. The matrix elements Fp ja in 
the row /3 can be substituted with Fp^ a ± aF 7i0 where 
7 is a row different from (3 and a is an element of Zk ■ 

— Row multiplication. Let a £ Zk- If a is not a factor of 
K, then all the matrix elements -F/3, a corresponding to 
the row (3 can be multiplied by a. 
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On a L x L torus, a basis for T{Q, K) is e m ,n(iu ii) = 
Si ltm 5i 2jTl . The homomorphism corresponding to the dis- 
crete Laplacian V 2 K can be represented as a L 2 x L 2 ma- 
trix. Before introducing the algorithm, we show that the 
system determining Ker(V^-) can be reduced by substi- 
tution from L 2 to 2L variables. The matrix form of the 
harmonic function h G Ti(Q, K) is 



/ 



h = 



/i(0,0) 
h(0,l) 
MO, 2) 



h(l,0) 

h(l,l) 

h(l,2) 



\h(0,L-l) h(l,L- 1) 



h(L- 1,0) \ 
h(L- 1,1) 
h(L - 1,2) 

/i(L - 1,L- l)/ 



(52) 



Equation ® yields h(j, 2) = 4/i(j, 1) - h(j - 1, 1) - + 
1, 1) — h(j, Q). In the same way the elements of the fourth 
row are determined by the elements of third and second 

row, and so on. This way the system V 2 K {h) = is reduced 
to a system of 2L unknown variables only, i.e. the elements 
of the two first rows of l|52() . This reduction procedure 
can be implemented into a computer program allowing 
to study the homomorphism kernels even for very large 
systems. 

D.2 Triangularization procedure 

First of all, we recall that Z p * is not a field, therefore 
we have to pay attention to some algebraic details. Any 
element a € Z p « can be decomposed as a = p h q with 
h < k and q an integer such that q mod p =/= (i.e. q 
admits the reciprocal q^ 1 S Z p » ), we will call this product 
p-factorization. Moreover, given two elements a and a' of 

Zp« with p-factorizations a = p h q and a' = p h q' , if b = a+ 
a' and c = p ha q c , we have that b = p hb qb with hi, > 
min(/i, /i'), and c = p' lc <7 c with h c = h' + hi (if (h + h') > k 

then c = 0). Taking into account these properties of Z pK , 
the matrix P 1 in equation (|50J) can be put in a triangular 
form a.s follows 

— Consider the p-factorization of all the matrix elements 
F/3, a — P hf> ' a< l0,a and find the elements Fp <a with the 
smallest exponent hp^ a . Swap rows and columns so 
that one of these elements is placed in the position 
(1,1) of the transformed matrix F A . Then Fp a = 

- Multiply the first row of Fp a by {qA^ 1 . We get 



P 



>ql a with qf fl = 1 and h? A < h% a V/3, a. 
Subtract to all the elements p hf3 - a q^ a (/? 7^ 1) the el- 



ements of the first row multiplied by p^c- 1 hl - 1 > 
i.e. transform F B into the matrix F% '„ = 



/3,a -P"^ a qp 



9/3,i; 
5 if 



/3 = landF|^/^ Q 
P^l.We have 



*.i- h i.i+*S«gfX a if 







V 







F 1 



(53) 



where F 1 is a matrix of size (M — 1) X (M — 1). As 
a consequence of the properties of sums and products 
of p-factorizations, Fl a — p h < 3 - a q\ a with hp > h\ A . 
The same property holds for Ff 2 , ■ ■ ■ 1 Ff M . 
— If all elements of F 1 equals zero the triangularization 
is concluded. Otherwise, apply to F 1 the same proce- 
dure applied to F. In this second case, we obtain a 
transformed matrix whose upper left element is given 
by p h2 ' 2 with h 2 ,2 > h A . 

The algorithm described above transform the general 
matrix F, into a triangular matrix of the form: 



/ T° 

^ 1 I 



F 



O 



TTTT 



7TTT 



E°- 



E°" 



E l 



(54) 



T K ) 



where T h (0 < h < re) are upper-triangular square matri- 
ces of size Q p h x Q p h, i.e. 



rph 



/ p h T( l . 

p h ■■■ T£ 



\ 



2,Q„h 



p h ) 



(55) 



therefore we have J2h=o Qp h = ^ where MxM is the size 
of the original matrix F. For T h we have Tp a — p mi3 - a qp^ a 
with Trig a > h. All the elements of the Q p h x Q p j rectan- 
gular matrices /tJ equal zero. While, for the Q v h x Q p j 



matrices E h ' j , Epf a 



P^qf*, with mf^ >fc. If 
ft = 1, if is a prime, A4(IC) is a vector space and not a 
module and we recover the usual Gauss elimination pro- 
cedure, in particular, all the elements have the reciprocal 
and we can discard the p-factorization. By substitution, 
one can directly check that, for each h, 1 < h < ft, equa- 
tion F^p ta k a = has Q p h independent solutions of pe- 
riodicity p h , which generate the submodule Ker(_F). 

The triangularization procedure for the discrete Lapla- 
cian can be directly implemented in a computer program. 
Main limitations consist in the computer capability in op- 
erations with integers and modules. For example, since one 
has to evaluate quantities such as (a-b) mod p K , p 2K has to 
be smaller than the maximum integer. Furthermore, the 
mod p K operation turns out to be very slow for large p K , 
limiting our calculations to values of p K up to 20000. On 
the other hand, the reduction from L 2 variable to 2L vari- 
ables, as described in the previous section, allows to study 
systems of a quite large size. As an instance, for a small 
value of p K (i.e. p — 103, k — 1), the kernel of V 2 K on a 
torus of size L = 400 can be calculated in a few seconds 
on a standard personal computer. 

The algorithm has been implemented in the follow- 
ing way. We considered all the primes p smaller than 
Kl = 20000 and we apply the algorithm to the opera- 



tor V 2 K for K 



P,P 



up to the value p K such that in 
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the diagonalization procedure we get only one function of 
periodicity p K i.e. the constant. In general we find that p K 
is smaller than Kl at least for not too large sizes. Theorem 
3 proves that it is useless considering larger p L . A simple 
analysis of the diagonal elements of the triangular matrix 
lf54"|) allows to evaluate how many elements of periodicity 
pi < pK. b e l on g S to Ig.Kg ■ By Lemma 2 such functions can 
then be written as elements oiH.(G,p L ). 

E Preflows and flows on graph 

Definition 4 Given an oriented graph Q , we call inte- 
ger "preflow" P(i.j) an integer-valued function {{i,j)} — * 
Z E , where Ne is the number of edges. The preflows con- 
stitutes a group V with the addition in Z. 

We expose here some properties of the integer preflow 
group V providing a interesting bridge between isoinvari- 
ant transformations and the graph topology. 

Definition 5 We denotes with {r} en the class of isoen- 
ergetic transformations t:Z^C^Z'^C defined as 
Z' = Z + t where t G Z n is an integer vector such that 
y^ - n = 0. Clearly {r} en forms an Abelian group. 

The isoinvariant transformations © {r/}i S o are a sub- 
group of {r} en . Moreover, we note that once it is fixed the 
state Z, Dcfinition[S]provides a one to one correspondence 
between the states of an energy surface and the elements 
of the group {r} en . In the same way the configurations of 
an atom are in a one to one correspondence with the isoin- 
variant transformations {rj}i S0 . Therefore, the number of 
atoms Af(G) can be evaluated as |{r} e „/{ry} iso | where | • | 
denotes the order of the quotient group. 

Definition 6 We denote with E the morphism associ- 
ating to each integer preflow Puj\ € V the isoenergetic 
transformation t = E(Puj\) : Z S C — > Z' € C , Z[ = Zi + 
Yj P{i,j) ■ We note that Range(S') = {r} en while we de- 
note as T = Ker(S') the subgroup of integer (conservative) 
flows, i.e. P( i;i) E T ifj^j p (i,j) = °- 

Theorem 5 S(Puj\) is isoinvariant if and only if Puj) 
belongs to the subgroup X of the irrotational elements of 
V , i.e. if for each loop of the graph io , ii , . . - , i n , io fik and 
ik+i are adjacent sites) P( io A 1 )+P(n,i 2 ) + - ■ -+ P (in,io) = °- 

Proof If S'(P(j ! j)) is isoinvariant we have • P(i,j) = (V 2 i7)i 
for a certain integer vector U. This means that Puj) = 
Ui—Uj yielding Pu,j) € T. Let P(i,j) € I, the integer vector 
U can be defined as U\ = and Ui — P(i,i t ) + P(u,i 2 ) + 
. . . + P(i ni i) (0, i\, . . . , i n , i is a path joining the vertices 
and i), the result does not depend on the choice of the 
path since P(i,j) is irrotational. We get J2j P(i,j) — (V 2 /7)i 
therefore S(P^ j)) is isoinvariant. □ 

We define as the dual space of V i.e. the set of all 
linear functions V — > Z. is isomorphic to V, indeed 
any linear functional can be written in a unique way as 
Yn «\ Pf \P(i 71 with P* £ V. 



Lemma 3 The dual space J 7 ^ is isomorphic to V/T. 

Proof The result follows showing that when Puj) € T 
then P$ d) Pm - if and only if P* }) el. D 

Theorem 6 The number of atoms generated by toppling 
invariants equals the graph complexity. 

Proof ^From Definition we get that {r} en is isomorphic 
to V IT. Since Theorem CD shows that X is isomorphic to 
Wiiso we get that M(Q) = \{T} en /{v}iso\ = \(P/T)/1\ = 
| (V /T) I T\ . Finally an important result proved in |17j shows 
that the graph complexity equals \T# / T\ and then from 
lemma|21 \ (J> jX) j T\ and this completes the proof. □ 
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